The Van der Waals interaction of the hydrogen molecule - 
an exact local energy density functional. 
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We verify that the van der Waals interaction and hence all dispersion interactions for the hydrogen 
molecule given by: 



W" = 



A_ 



B_ 



C 



(1) 



in which R is the internuclear separation, are exactly soluble. The constants A = 6.4990267..., 
B = 124.3990835 . . . and C = 1135.2140398 ... (in Hartree units) first obtained approximately by 
Pauling and Beach (PB) [Q] using a linear variational method, can be shown to be obtainable to 
any desired accuracy via our exact solution. In addition we shall show that a local energy density 
functional can be obtained, whose variational solution rederives the exact solution for this problem. 
This demonstrates explicitly that a static local density functional theory exists for this system. 
We conclude with remarks about generalising the method to other hydrogenic systems and also to 
helium. 

results. We shall show that the method is generalisable 
to other hydrogenic systems and also to helium for which 
systematic approximations can be derived. In section ||, 
we shall discuss the background to the exact solution ver- 
ifying that the PB variational method is essentially exact, 
though more slowly convergent than the SK method. In 
section [1] we shall present our method for solving the 
SK equation and show our results. In section |^ we shall 
formulate the local density functional theory (LDFT), 
discuss its solution and show that it converges to the ex- 
act results of section [11 . In section ^ we shall discuss 
the problem of helium and conclude in section VI with 
discussions about further work. 



I. INTRODUCTION 

Amongst the few non-trivial many-body problems in 
quantum mechanics, the hydrogen molecule was the first 
system to be thoroughly studied and continues to be re- 
searched, motivated by experiments ||^, Q and by 
advances in modern LDA techniques of computation . 
It is perhaps not so well known that the dispersion forces 
for this elementary system is amenable to an exact so- 
lution, due to certain confusion in the early literature 
with regard to the methods used to attack it. In this 
paper we shall use a method first propounded by Slater 
and Kirkwood (SK) |^], but whose equation (see eqn(||) 
below) they were unable to solve, to show that the dis- 
persion forces for the hydrogen molecule are exactly solu- 
ble. This method later became formalised as the method 
of Dalgarno and Lewis and is particularly suited for 
the problem of dispersion forces in general. However it 
does not seem to have been exploited in recent studies 
of the van der Waals interaction using density functional 
theories, 1^ . We shall demonstrate the superior con- 
vergence properties of the SK method. This was already 
known to the early pioneers, in contrast to fre- 

quency dependent methods, essentially based on a sum- 
mation over dipole matrix elements with excited states or 
an integration over the dynamical polarizabilities. The 
latter was derived through the original work of Eisen- 
chitz and Londo n [[10| , and had a rather strong influence 
on later studies ||11[] . For the hydrogen system, we shall 
derive an exact local density functional theory, soluble 
for this case, and whose solution converges to the exact 



II. BACKGROUND 

The Hamiltonian for the dispersion forces of the hy- 
drogen molecule was first derived by Margenau ||l2| and 
has appeared in subsequent editions of many textbooks, 
especially the famous one of Pauling and Wilson [p^ . 
The latter contains an excellent survey of the variational 
treatment for the van der Waals interaction for the hydro- 
gen molecule and also to early results for the helium sys- 
tem, which to this author's knowledge has not yet been 
updated. As is now well known, the form of the Hamil- 
tonian, derives essentially from a large distance expan- 
sion of the electron-electron interaction for the hydrogen 
molecule. This is given by | |T^ : 

H' = {l/R'^){xiX2 + 2/12/2 - 2ziZ2) + 

+ ^il/R^)[rlz,~rlz, 

+ {2xiX2 + 22/12/2 - 'iziZ2){zi - Z2)] 

A{^lR'')\r\Tl-^rlz\-hT\zl 



15z-i^Z2 + 2{xiX2 + J/12/2 + 42122) ] 



(2) 



In this expression xi,j/i,zi are the Cartesian coordi- 
nates of electron 1 relative to its nucleus, while X2,y2, ^2 
are those of electron 2 relative to its own nucleus 
and the z axis is directed from one nucleus to the 
other. The expressions contain the dipole-dipole inter- 
action (first term), the dipole-quadrupole interaction 
(second term) and the quadrupole-quadrupole inter- 
action (third term), with the first term being the well 
known van der Waals attraction which is dominant for 
large distances R - the internuclear separation. Through- 
out this paper we shall be using Hartrec units in which 
{e — h = ao — 1) and thus the absolute ground state en- 
ergy Eq which we shall require later has the unit of 1/2 
Hartree. In view of the symmetry of the various terms, 
it can shown [Q that the above eqn(H) is equivalent with 
respect to its second-order perturbation energy to the 
Hamiltonian: 



i7' = -2[a ^16 cos 6*1 cos 6*2 
+ /3 Ci?|cos6'i(3cos^ 6*2 - 1) 
+ 7 C??2 (3 cos^ 01 - 1 ) (3 cos^ 62 



1) 



(3) 



in polar coordinates, whereby: ^1.2 — 2ri^2,ct — 

{6)i/8 = (30)*/32 R^'^ and 7 = (70)*/128 R'^. 

In this section we shall briefiy review some of the difficul- 
ties in an accurate determination of the constants A, B, 
and C. One standard formula is the direct summation 
method as originally used by Eisenchitz and London [0 , 
for the van der Waals energy : 



12 



fl,nfl 



(i-;r^)(i 



r)(2- 



(4) 



where fi^m — 2(£'m — Ei)\zi,m\'^ is the oscillator strength 
as defined in terms of the dipole matrix elements zi_rn 
between the states {I, m) However, the series eqn(^) 
converges badly even for the discrete states and there are 
terms involving the matrix elements between discrete/ 
continuum and continuum/continuum states that are dif- 
ficult to evaluate but which ultimately determine the ac- 
curacy of the result. The value for the constant A « 6.47 
originally given by Eisenchitz and London ||lotl after much 
work testifies to the difficulty of this approach. Never- 
theless, eqn(||) has its appeal in that it can be recast into 
the form of an integral over imaginary frequencies of the 
dynamical polarizabilities of the two atoms: 



o roc 



OLi{iCia2{iOi 



(5) 



from which current theories for the van der Waals in- 
teraction for more complex systems like He are based. 
These employ a combination of linear response and time- 
dependent density functional theories to obtain the polar- 
izabilities . An integral equation is ultimately involved 
in the .solution for (Ti ffenerallv involving' various 



decoupling approximations, and thereafter the integral 
eqn(H) has to be performed. To the best of this author's 
knowledge, none of the theories proposed so far have been 
tested with the exactly soluble case of the H2 molecule 
Q , 1^ providing added motivation for our present work. 

Let us mention at the outset that the central difficulty 
of this problem has to do with an accurate treatment 
of excited states, their matrix elements with the ground 
state (both discrete and continuous) and the poor con- 
vergence of the series like eqn(^). These are all sticky 
points with the modern local density functional theo- 
ries (LDFT) 0. It has been identified by Slater and 
Kirkwood (SK) in their classic paper j^, that a supe- 
rior method involves a direct perturbation wavefunction 
ansatz of the form: 



'0(ri,r2) ■(/'o(r-i,r2)[H-0(ri,r2)], 



(6) 



in which is the ground state wavefunction of the un- 
perturbed system. The function (a two particle correla- 
tion function as we shall see), satisfies the following exact 
differential equation, easily derived from the Schrodinger 
equation up to first order in v which can be va,vb or 
vc accordingly. Following the notation of SK Q , this is 
given by: 



1 



V^(/) -I- (VlnV'o)-(V(/)) -v = Q. 



(7) 



SK used this equation as their basis for the treatment of 
the H2 and He systems. For the moment we shall con- 
centrate on the H2 molecule, which by the substitution: 



(8) 



leads to the differential equation (hereafter known as the 
SK equation) in the case of the van der Waals interaction 
VA as given by: 



d^R d^R 



+ (7-1) 



R{ 



dR 4 



1 



dR 



0. 



(9) 



Note that i?(^, ^') is strictly non-local, but it can in prin- 
ciple be derived from a local density, as we shall see. Un- 
fortunately SK were unable to solve this inhomogeneous 
(PDE) equation and resorted to various approximations, 
one of which is to ignore the differential terms and as- 
sume: 



4C + e 



(10) 



They suggested that this is a good first approximation 
and that subsequent approximations can be obtained by 
substituting this into the differential function of eqndfl) 
and iterating. They evaluated the constant A « 6.14 



using eqn(px|), a result tabulated in the book of Paul- 
ing and Wilson fl^ , but unfortunately we have found 
this value to be in error. The correct approximate value 
should be A « 6.23. The second error in the SK paper 
is that their proposed iteration method does not work. 
In fact it is a non-convergent procedure, and attempts 
to employ it leads to divergent results. Earlier on in 
this work the author has carried out an iteration of their 
scheme to four orders and found the results diverging. 
Nevertheless SK suggested other approximation schemes 
like the ansatz Ar^r'" where A and v are variational pa- 
rameters (see later) which they have found to give good 
results for both the H2 and He systems. We shall discuss 
the solution of eqn(^ later in the next section. Here 
we shall mention the best solution for the H2 problem 
to date. This was the widely cited paper of Pauling 
and Beach (PB) Their solution employs a general 
variational method in which the matrix elements for the 
Hamiltonian were evaluated using special orthogonal or- 
bitals constructed for a solution to the Stark effect prob- 
lem They have found these orbitals to be ideal for 
this problem by which all matrix elements for the inter- 
actions VA,B.c can be computed accordingly. Thereupon 
they were able to set up an infinite determinant from the 
secular equation which they have evaluated up to a rank 
of order (26 x 26) to find the energies. Their results for 
A = 6.49903, B = 124.399 and C = 1135.21 (accurate 
up to the last decimal) remain the most accurate to date, 
until our present work and is very impressive for 1935. 
However they cautioned that their "treatment has not 
led to an exact solution" owing to the uncertain nature 
of their variational method and the convergence proper- 
ties of their wavefunctions fl^ ] . In addition their method 
does not yield the expansion coefficients for the wave- 
functions, needed to ensure normalisability and thereby 
obtain the normalisation constant C{R). Furthermore 
their procedure cannot obtain the perturbation in the 
charge density which will be of interest to us here. The 
reader is referred to their earlier paper for information [Q , 
many of whose details are now merely of historical inter- 
ests. Neverthless, they have identified their method as 
identical with and is a more generalised variational form 
of that used by Hasse ||l^ , who also was the first to treat 
the H2 and He problems. In the next section we shall dis- 
cuss the exact solution that SK had failed to obtain for 
eqn(^). We should note that the power of their method 
lies in the ability to employ the interaction function it- 
self to project out all relevant components of the excited 
states for second order perturbation theory, into the two 
particle correlation function R{r,r'). This then has a 
concise form satisfying an inhomogeneous PDE eqn(||). 
This method was later generalised to higher order per- 
turbation theory by Dalgarno and Lewis |^] and Schwartz 
sec also Schiff Wh- 



III. SOLUTION OF THE SK EQUATION 

Having set the background to the hydrogen problem we 
shall next discuss the solution of the SK equation eqn(H). 
This is a two dimensional inhomogeneous (PDE) and we 
must first start by discussing the appropriate boundary 
conditions. This boundary value problem is unusual in 
that it is not of the standard Dirichlet or Neumann type 
as is common in electrostatics In fact there are 

no particular a prior boundary values apart from the 
requirement for the normalisability of the wavefunction, 
and special values such as i?(0, 0) are determined only 
after the solution is obtained. Nevertheless an analogous 
Green's function integral equation method which is exact 
is known to exist due to the work of Levi ||l9|. This is 
of no interest to us here, so that we shall merely outline 
the method in an appendix, but it may be useful for 
making connections with other integral equation methods 
of treating the problem, such as via linear response theory 

i, §■ 

The method we shall use to solve eqn(^) is gained from 
experience in solving the two sphere problem of classical 
electrostatics, |20|, |2l|]. As we shall see the convergence 
of this problem by our method is superior to the elec- 
trostatic case of two spheres, which required nearly 200 
terms for convergence to only two decimal places . We 
begin by expanding the two particle correlation function 
in the following ansatz in terms of orthogonal polynomi- 
als: 



(11) 



where the functions Hn(x) = L^^j^(a;) are defined in 
terms of associated Laguerre polynomials . Inserting 
this into the SK eqn(||), it reduces to the form: 



E 



a-L: 



(l-l) , (n-1) 



1 

4' 



(12) 



when use is made of the equation for Laguerre polyno- 
mials. Upon multiplying both sides of eqn([l^) by: 



(13) 



and then integrating, with the use of the properties of the 
Laguerre polynomials, we derive the following infinite set 
of linear equations for the am^s- 

- am^si^sgrnQs + 2mgsqm) + am,s-i(s + '^)gmqs-i 

(.-1) 

Hm — 1 

+ a„i+i^s(m - l)gsqm+i = ^m,s- (14) 
where: 

1 (s-l)(s + l)! (s + 1)! 

9s = 777 1 TK-f ' 'i'^ = 



144 



(,s-2)! 



(.-2)!' 



("15^ 



This set of equations is readily solved symbolically using 
Mathematica version 3.1 on a PC. We shall present the 
results in the next subsection. Here we shall obtain the 
form of the energy expressions. Firstly, as was shown 
first by SK for an arbitrary correlation function R 
not necessarily satisfying eqn(^, the energy expression 
is given by: 



1 



(16) 



where C is the differential operator given by the LHS 
of eqn(^). This form is particuarly useful when we 
look at density functional theories afterwards. At this 
point we shall merely mention that the change neces- 
sary for calculating the dipole-quadrupole energy es and 
the quadrupole-quadrupole energy ec are the form of v 
which becomes vb and vc respectively. The form of R 
also changes which we shall call Rb and Re, (subscripts 
being only used when there is a need to avoid confusion) 
and so are the operators Cb and Cc ■ They can be ob- 
tained from the appropriate SK equations and are given 
by: 



d^RB 



d^RB A .BRb ,6 .ORb 

-i?B(i + |)- J = 0, (17) 



and 



d^Rc , d-'Rc 



.6 dRc 



.6 dRc 



2i?c(^ + ^)- J = 0. (18) 



They can be solved by the same method as for Ra, the 
expressions for the expansion eqn(|ll[) now being: 



(19) 



RB{^,e) = J2^i,nHi{0Gnie) 



l,n 



and 



(20) 



where G„(a;) = L^_^_2{x) is a higher order associated La- 
guerre polynomial. We note that for the exact solution 
C[R{^,£^')] — from which the energy expressions are 
easily obtained in terms of the first few expansion co- 
efficients. We shall collect the formulas for the various 
energy constants in terms of these coefficients as: 



B = 

n 



-12[a2,2 - 02,3 - ^3,2 + 03,3], 

-270[62,3 - hA ~ &3,3 + ^3.4], 
= — 2835 fr.Q 5 — r.Q a — c.a 1 4- r.A a}. 



Thus the energy constants can be determined to any de- 
sired accuracy using symbolic manipulation codes such 
as by Mathematica version 3.1 which yields the a,b,c 
coefficients as exact fractions. 

In order to check the definite convergence of the wave- 
functions, a point of concern for Pauling and Beach [Q, 
we have also computed the normalization constants for 
the wavefunctions in eqn(^. These are given by the in- 
tegrals: 



a{R) 



^0 



dr / dY'%vtR. 



)2 



(22) 



where i = A, B,C respectively. More appropriate ex- 
pressions are given in terms of the constants Di where 
we have factored out the distance dependence R: 



D 



(23) 



The values of Di can thus be obtained in a series in terms 
of the a,b,c coefficients and computed to any desired ac- 
curacy. Of particular interest to us is a calculation of 
the density. This can be obtained by direct partial in- 
tegration of the wavefunction in eqn(^. Since we are 
only interested in the density perturbation, by subtract- 
ing the unperturbed density po — ne^^ for the ground 
state hydrogen atom, we have: 

Me,e) p(e,^)-po 



Po 



Po 



^6*' sin 61' / d^e-^'C^ 
Jo 

v\^,^',0,d')R\ta- 



(24) 



We shall redefine alternative functions /i(C) which we 
shall call "densities" and they are the main focus in this 
paper: 



SpA{^,e) 

Po 



Po 

5p%\^,0) 
Po 

Spcite) 

Po 



e cos' Of 



16i?8 



^ e^(3cos^^-l)Vf (0 



64i?8 
7 



256i? 



10 



(25) 



Note that there are two densities for B, since the dipole- 
quadrupole interaction is asymmetric. All these densities 
are readily computed from the appropriate wavefunction 
coefficients. For example we have: 



MO = ^n,mHn{OHm{£,) CtC, (26) 

where: 

I 

- [I - l)g/+ia„i,/+i], (27) 

and so on. Their resuhs will be given in the next subsec- 
tion. 



A. Exact results 

The infinite set of linear equations such as eqn(p^ is 
truncated at each order and the solution for the coef- 
ficients a, 6, c are solved by the use of Mathematica on 
a PC accordingly. It is remarkable that the results are 
so fast converging unlike that for the two spheres prob- 
lem in electrostatics . Previously for the two spheres 
we have found the need to export the codes to a Silicon 
Graphics workstation running Mathematica, as expan- 
sions up to the order of 200 coefficients are necessary 
before we could obtain convergence to 2 decimal places. 
The present codes run readily on a PC and as the follow- 
ing Table || shows they converge rapidly. The superiority 
of convergence of our approach versus other methods |lC|] 
is obvious from the results shown in this table. Our con- 
vergence rate is even better than the calculations of PB 

0- 

TABLE I. Computed constants A and Da 
order of truncation A Da 



1 


6 


6 


2 


6.22222222. . . 


6.61728395. 


3 


6.46153846. . . 


6.88757396. 


4 


6.48214285. . . 


7.40242346. 


5 


6.49844398. . . 


7.40024688. 


6 


6.49900257. . . 


7.39872679. 


7 


6.49902535. . . 


7.39863094. 


8 


6.49902659. . . 


7.39862559. 


9 


6.49902669. . . 


7.39862525. 


10 


6.49902670. . . 


7.39862522. 



TABLE II. Computed constants B and Db 
order of truncation B Db 



1 


115.71428571. . . 


24 


79591836. . 


2 


118.96875 


26 


74423828. . 


3 


124.26672692. . . 


30 


14754412. . 


4 


124.39502505. . . 


30 


12987113. . 


5 


124.39891831. . . 


30 


12699933. . 


6 


124.39907397. . . 


30 


12683881. . 


7 


124.39908277. . . 


30 


12682990. . 


8 


124.39908349. . . 


30 


12682930. . 


9 


124.39908357. . . 


30 


12682924. . 


10 


124.39908358. . . 


30 


12682924. . 



Here we see that at truncation order 10, where there 
are about 100 coefheients, the smahest being aio.io, 
we have achieved convergence to at least seven decimal 
places. Our results are in exact agreement with Pauling 
and Beach |jl[ , indicating that they have indeed found the 
exact energy variationally. In particular the normalisa- 
tion constants converge to the same accuracy indicating 
that the wavefunctions are normalisable and thus well be- 
haved. In the following Figj^ we shall show the computed 
density function 



f 




FIG. 1. Computed density function for successive 

truncations from order 5 (dashes). Note the convergence in 
the successive plots towards the solid curve. 



TABLE III. Computed constants C and Dc 



order of truncation 


C 


Dc 


1 


1063.125 


199.3359375 


2 


1132.610294117. . . 


238.583207179. . . 


3 


1135.107421875 


238.686733245. . . 


4 


1135.208820466. . . 


238.645331783. . . 


5 


1135.213725627. . . 


238.641796376. . . 


6 


1135.214015982. . . 


238.641543683. . . 


7 


1135.214037581. . . 


238.641524775. . . 


8 


1135.214039617. . . 


238.641523160. . . 


9 


1135.214039858. . . 


238.641522996. . . 


10 


1135.214039892. . . 


238.641522976. . . 



It is to be noted that the large distance behaviour in 
the density dictates the subsequent accuracy in the en- 
ergy calculation. We show the density at various orders 
in Fig.| 



f 




5 10 15 20 25 30 

FIG. 2. Computed density function a-t larger dis- 

tances for successive truncations. Order 5 is the bottom curve 
(dashes) while order six is the top curve (line). The curves 
have alternative curvatures for odd and even orders of trun- 
cation. Note the convergence to the middle plots at higher 
orders and the difference with the scale of Fig.|^. 

In a similar way we have computed the B and C energy 
constants in eqn(gj), as well as their respective normalisa- 
tion constants Db and Dc respectively. These are shown 
in Tables || and |ll[ Note that the numbers are tabulated 
as decimals for convenience, but they are exact fractions 
as given by Mathematica. As such some of the numbers 
terminate as decimals, whereas the others have a finite 
period, most of which are much longer than the 8 or 9 
decimals shown. For completeness of the results we have 
plotted all the various density functions calculated from 
the exact wavefunction coefRcients. 
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FIG. 3. Computed density function f^g\Q for successive 
truncations from order 5 (dashes), as in Fig.nl 



FIG. 4. Computed density function fg\S,) at larger dis- 
tances, for successive truncations from order 5 (bottom 
curve), as in Figj^. 



FIG. 7. Computed density function /c(C) for successive 
truncations from order 5 (dashes), as in Fig.pl . 




FIG. 5. Computed density function fg^Q for successive 
truncations from order 5 (dashes), as in Fig.l 




FIG. 6. Computed density function fg {^) at larger dis- 
tances for successive truncations from order 5 (bottom curve), 
as in Fig.W 




5 10 15 20 25 30 

FIG. 8. Computed density function /c(C) larger dis- 
tances for successive truncations from order 5 (bottom curve), 
as in Fig.^ 

Note that for shorter distances in figures |^J|,|^ and 
0, the convergence proceeds monotonically from the top 
curve to the bottom. These plots are nearly straight, but 
for larger distances they have alternative curvatures for 
odd and even truncation orders. These observations fur- 
nish interesting approximation methods, which will not 
be investigated here. In the next section, we shall use 
the results obtained so far to derive an exact local den- 
sity functional theory for this system. 



IV. AN EXACT ENERGY LOCAL DENSITY 
FUNCTIONAL 

Current thinking in density functional theory (DFT) 
approaches the problem of dispersion forces from eqn(||), 
as mentioned earlier, via time-dependent generalisations 
of DFT. It is interesting to note that the method of 
Kohn et al Q has yielded results in agreement with the 
best theoretical value for He up to two decimal places, 
in spite of a 7% error in the completeness sum rule in 
their calculations. However it is difficult from these the- 
ories to develop systematic improvement methods and 
would require considerable expertise in time-dependent 



density functional methods. The question arises a long 
ago in a paper pointed out by Lieb |23|| , that the universal 
Hohenberg-Kohn (HK) functional I[p] has "hidden 
complexities" with respect to the van der Waals interac- 
tion. The dynamical dipole fluctuation properties leading 
to the latter "is somehow built into I[p\, but an explicit 
form of I[p] that will produce this effect has yet to be 
displayed." This is in spite of a very general proof for the 
universal character of van der Waals forces for Coulomb 
systems [ psf . We shall answer this question some way 
for the H2 molecule. It is noteworthy that perhaps as a 
result of these remarks, the search for an accurate static 
local density functional theory for dispersion forces has 
more or less been abandoned. In this section we shall 
demonstrate that for the H2 molecule system at hand we 
can formulate an exact static DFT. We start from the 
observation that the two particle correlation R(r,r') is a 
functional of the density Ja- In eqn(p6|) «/the quantities 
ctn,m are fixed, thereby fixing the density, then eqn(p7|) in 
principle can be inverted to obtain the coefficients a„^m 
thereby determining R[r,r'). Therefore upon substitut- 
ing this R[an,m] into the energy functional eqn([l6|), then 
a variation of ca with respect to ari,m would yield the 
exact ground state energy. This corresponds to the con- 
straint search algorithm of Levy , but note that this 
is not the HK functional as it is specific to this problem. 
Formally we can write this as: 

{mTTEo) j„ j„ 

Then the ground state energy and density can be ob- 
tained from: 



SeAifAim 



0. 



(29) 



However eqn(g^ is not the most convenient to use. Its 
inversion corresponds to a non-linear programming prob- 
lem with many solutions. This complexity came from 
our choice of orbitals which for the exact solution fixes: 
£[i?[/yi(^)]] = 0. An alternative choice of orbitals can be 
made which then fixes ^[/a(C)] but now £[.R[/a(^)]] 0. 
Nevertheless eqn( p9h w ill still yield the exact ground 
state energy via eqn(E8h which is the essence of our DFT. 

A convenient choice of orbitals is determined from the 
density expression: 



fAiO 



d^'e-«'ri?i(e,e')- 



(30) 



It can be easily seen that the choice of orbitals (f>n{0 
Ll{£) such that: 



iOMO etc. 



where: 



an,i 



m 



i'^«(C)0m(C'): 



(32) 



(33) 



in which LUm = m\/(rn — 4)!, as appropriate for these 
orbitals. The variation of the density /a can now be 
effected by directly varying the coefficients dn,m- The 
latter can be easily computed from eqn(p^), which can 
be carried out symbolically as well. The integrals re- 
quired throughout the calculation can also be computed 
in closed form, facilitated by the symbolic integration 
capabilities of Mathematica. Our results are tabulated 
in Table [V. The results are in exact agreement with 
Table |, since as we have noted earlier the output are 
exact numeric fractions that can be co mpare d order by 



order with the results of the subsection HI A . The inte- 



grals involved are somewhat more lengthy here so that 
we have not extended the calculations beyond the 6th or- 
der. In the following figures, we have plotted the density 
functions which are again in exact agreement with the 
previous figures |^ and ^. 
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FIG. 9. Computed density function /a(0 truncated at or- 
der six, the dotted curve is our exact DFT result and the solid 
curve is the exact result of the previous section. Both curves 
are identical. 



TABLE IV. Computed constant A using our DFT eqnfcg 



order of truncation 



A 



6.46153846. 
6.48214285. 
6.49844398. 
6.49900257. 



(31) 



nrovide a much simnler form for the derisitv whereunon: 
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FIG. 10. Computed density function for truncations 

at order 6 for larger distances. The dash curve is our ex- 
act DFT result and the solid curve is the exact result of the 
previous section, both curves are seen to be identical. 

We can show in the same way using the exact DFT 
for the other energy constants B an d C that they yield 
identical results as subsection [II A. The convenient or- 
bitals for these other density expansions are easily seen 
to be: 0„(OXm(e) for fi'\0 and x« (0x^(0 for '(C) 
and /p^^(^) accordingly, where Xn(C) = -^n(C)- We found 
these results to be very instructive from the viewpoint of 
DFT. In particular approximate densities can be devel- 
oped. For example the form: 



/^(O- Const e 



2iy 



(34) 



with Const = A (4 + 2;/)! follows from the SK ansatz 
for R{r,r') = Xr'^r"^. The variational results using this 
approximation (which have to be computed numerically) 
give two significant figures accuracy except for the case 
of es, see Table 0. 

TABLE V. Computed constants A,B and C using eqn(Q 



A 


B 


C 


6.48965 


116.795 


1134.71 



The following figures compare the approximate with 
the exact densities. We shall not discuss these approxi- 
mations here as detail investigations will require further 
work. In the next section we shall discuss if our results 
could be extended as approximate methods for more com- 
plex systems for which no exact solutions are known. 
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FIG. 11. Approximate density function /a (C) (dashed) cal- 
culated from eqn(p4|) compared with the exact one. 
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FIG. 12. Approximate density function /s(5) (dashed) cal- 
culated from eqn(|34[) compared with the exact one. Note that 
this poor density also yields a poor energy constant. 
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FIG. 13. Approximate density function fc{£,) (dashed) cal- 
culated from eqn(|34[) compared with the exact one. Note that 
this good density also yields a good energy constant. 



V. HYDROGENIC SYSTEMS AND HELIUM 

For hydrogenic atoms for which we can replace the core 
charge from Z — 1 by Z' say, the modifications of our 
method is quite straightforward. A first principles calcu- 
lation however is a different matter. We shall only discuss 
the van der Waals energy ca from hereon. This can be 
seen by considering the helium problem. The interac- 
tion energy contains a sum of interactions taken from all 
possible pairs of electrons between the two atoms: 



N N 
j = l k=l 



(35) 



where N is the number of electrons on each atom and: 



(36) 



The two particle correlation function eqn(^ now breaks 
up into a sum of pairs: 



N N 



'/'(I'j' i-fc) ^ ^j-kR{rj,rk), (37) 

" j=i k=i 

and the SK eqn(|) now becomes a set of N"^ equations: 



d^R d^R ,4 d\nija,dR ,4 dlni^Q.dR 



d^k d£,k 
(38) 



These equations are now coupled, since in general ■00 
is a many-body wavefunction. Hence the problem is in 
general insoluble. Nevertheless, the situation in which 
Tpo is given in an LDA approximation as a product of 
Kohn-Sham (KS) orbitals |^ is greatly simplified and 
should form the basis of an LDA approach as we shall 
see . We have learnt from the hydrogen problem that 
the correlation function R{rj,rk) and the density f{rj) 
are closely connected with these orbitals. For hydrogenic 
atoms in which the core can be considered as a closed 
shell, then the following approximation: 



9 In -00 



-Z'/2, 



(39) 



where Z' < 1 can be made without a significant lost of ac- 
curacy. In this case eqn(p8|) reduces to a single equation 
and is amenable to analytical treatment as for hydrogen. 
However from the form of the energy, which is additive 
in terms: 



TV TV 
j=l fe=l 



(40) 



1 I ^lkR,,ka - ^.ARj.kMdr 



(41) 



we will be motivated to consider a DFT theory such as 
that presented in section IV with appropriate approxi- 
mations. Note that the integral in eqn(^l|) is over the 
coordinates of all the electrons with the implicit depen- 
dence of Rj,k on the others. As can be easily seen, the 
operator Cj^k simplifies considerably if the ground state 
of the atom tpg is well approximated by a Hartree type, 
or KS type wavefunction for spherical atoms: 



^9 



TV 



(42) 



as in this case the eqns(^8|) decouple. The accuracy of 
the calculations will be dependent on approximations to 
the density f{rj) and the wavefunction ^q. Further in- 
vestigations along these lines will be able to provide a 
systematic study of van der Waals interactions as in the 
case of hydrogen detailed here. 



VI. CONCLUSION 

This paper sets out an exact solution for the van 
der Waals and other dispersion forces for the hydrogen 
molecule using the method of Slater and Kirkwood, 
By considering the density distributions fA,B,c we have 
shown that in this case, an exact energy density func- 
tional exists for this problem which when minimised with 
respect to the density, yields the exact results. We have 
shown that the energy constants A, B, C can be calcu- 
lated to any desired accuracy, the first few decimals being 
in full agreement with Pauling and Beach We have 
also considered the extension of our method for more 
complex systems such as hydrogenic systems and helium 
for which approximations must be invoked. A systematic 
study of these and generalizations to include the effect of 
a surface will be the subject of future work. 

Acknowledgements: The author wish to thank the 
NCTS (Taiwan) for their hospitality during the final 
stages of this work. This paper is dedicated to the mem- 
ory of J. Mahanty, who crossed path with the author 
during a period at the Australian National University in 
the late 1980's. 

Appendix 
Integral equation method 

The integral equation technique for solving the SK 
eqn(^) is due to Levi ||l^. We define the operator W 
acting on any function such as R: 



W[R] = AR + aR^ + bR^> + cR, 



(43) 



where: 



where the subscripts denote differentiation from hereon 
and A is the two dimensional Laplacian operator, so that 
the SK eqn(^ is given by: 



in which: 



W[R] = /; 



(44) 



: 1. 



(45) 



A uniqueness theorem for particular solutions can be 
proved for any value of c < 0, thus guaranteeing 
a solution for eqn(^^. With the use of an appropriate 
two dimensional Green's function G(x, y\x' , y') such that: 



AGix, y\x', y') = -2n5{x - x')5{y - y'), 



(46) 



then it can be easily shown that for any arbitrary function 
uj{x, y) the solution R{x, y) is given by: 

R{x,y) = Lo{x,y) 

+ fdx' [ dy'G{x,y\x',y')p{x',y'). (47) 



The function p{x, y) is given by the solution of the inte- 
gral equation: 



P{x^v) j dx' j dyK{x,y\x\y')p{x',y')+g{x,y), 

(48) 

where the kernel K{x, y\x' , y') is of the form: 



K{x,y\x',y')^^(^aG. 
and the function 17(2;, y) is: 



bGy + cG 



(49) 



(50) 



That eqns(47) to ( pOD give a solution for eqn(44) can 
be easily shown by operating on R in eqn(^) with the 
operator W. With an appropriate choice of Green's func- 
tion G{x,y\x' ,y') and uj{x,y), the iteration of eqn(|4^) is 
equivalent to our solution for R{x, y) as obtained in sec- 
tion fll. 
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